%% model_calib.m
%
% Set up model calibration
%--------------------------------------------------------------------------


%--------------------------------------------------------------------------
%% Discount Function Calibration

s = 2; % risk aversion in CRRA utility
switch calibrationCase
    case 1
        %Exponential Benchmark (No Slowness)
        beta = 1;
        rho = 0.0125;
    case 2
        %Present Bias Benchmark (Procrastination)
        beta = 0.83;
        rho = 0.0088;
    case 3
        %Exponential (Procrastination)
        beta = 0.999;
        rho = 0.0133; 
    case 4
        %Exponential (Rationally Slow)
        beta = 1;
        rho = 0.0134;
        slow_refi = 1;
        rational_slow_refi = 1;
    case 5
        %Present Bias (No Slowness)
        beta = 0.675;
        rho = 0.0054;
        slow_refi = 0;
    case 6
        %Present Bias (Rationally Slow)
        beta = 0.725;
        rho = 0.007;
        slow_refi = 1;
        rational_slow_refi = 1;
end


%--------------------------------------------------------------------------
%% Housing
if ~exist('h', 'var')       
    h = 3.29;               %house size
end
thetam = 0.8;               %max LTV
nocashoutFlag = 0;          %Dont allow for cash-out refis (robustness)
ARMFlag = 0;                %ARMS instead of FRMS          (robustness)

ka_rewrite = 0.05;          %fixed cost of adjustment
ka_prepay = 0.002;
    assert(ka_rewrite >= ka_prepay);
mortgagehalflife = 20;
    xi = -log(0.5)/mortgagehalflife;

%Refi Slowness
la_slowrefi = -log(0.5);
    if ~exist('slow_refi', 'var')
        slow_refi = 0;
        if beta < 1
            slow_refi = 1;
        end
    end
    if ~exist('rational_slow_refi', 'var')
        rational_slow_refi = 0;
    end
    if rational_slow_refi == 1
        assert(slow_refi == 1, 'Must set slow_refi = 1 in this case');
    end


%--------------------------------------------------------------------------    
%% Liquid asset state space
Nb    = 251;                %number of grid points
bmin    = -1/3;
bmax    = 2.5;
if bmin == 0
    b = linspace(bmin,bmax,Nb);
    db = b(2)-b(1);
    bzeroInd = 1;
elseif bmin < 0
    %ensure that 0 included in b-grid
    ratio = abs(bmin)/(bmax-bmin);
    lowerPts = ceil(Nb*ratio);
    b = linspace(bmin,0,lowerPts+1);
        db = b(2)-b(1);
    b = [b, [db:db:db*(Nb-lowerPts-1)]];
    bmax = max(b);
    bzeroInd = find(b == 0);
end


%--------------------------------------------------------------------------
%% Mortgage state space
Na   = 61;                  %number of grid points
amin   = 0;
amax   = thetam*h;
a  = linspace(amin,amax,Na);
da = a(2)-a(1);

%Adjustment decision set-up
NaP_rewrite = 200;
NaP_prepay = 100;


%--------------------------------------------------------------------------
%% Income Process
averageIncomeSCF = 95718;
Nz = 3;

%Guerrieri-Lorenzoni Calibration (+ Floden and Linde)
Delta = 1/4;
rate_y_OU = -log(0.967)/Delta;
sig_y_OU = sqrt((0.017*2*rate_y_OU)/(1-exp(-2*rate_y_OU*Delta)));

%Construct transition matrix (imposing reflecting barrier)
dz = sig_y_OU;
la_mat_z = [-rate_y_OU*dz/dz - (sig_y_OU^2)/(dz^2) + (sig_y_OU^2)/(2*dz^2), rate_y_OU*dz/dz + (sig_y_OU^2)/(2*dz^2), 0; ...
            (sig_y_OU^2)/(2*dz^2), -(sig_y_OU^2)/(dz^2), (sig_y_OU^2)/(2*dz^2);...
            0, rate_y_OU*dz/dz + (sig_y_OU^2)/(2*dz^2), -rate_y_OU*dz/dz - (sig_y_OU^2)/(dz^2) + (sig_y_OU^2)/(2*dz^2)];
assert(min(round(sum(la_mat_z,2), 14) == zeros(Nz,1)) == 1)

%Income states
zstates = exp([-sig_y_OU, 0, sig_y_OU]);
    %Adjust so mean income = 1
    [weights,~] = eigs(la_mat_z',1,'lr');
    zstates = zstates.*(sum(weights)/(zstates*weights));
    
%Fiscal policy shock size ($1000, scaled by SCF income)
liqShock = 1000/averageIncomeSCF;
mortgageShock = -liqShock;


%--------------------------------------------------------------------------
%% Interest Rate Process
Nr = 4;
rbAll = [-0.01, 0, 0.01, 0.02];

rmean = mean(rbAll);
rbdeath = rmean;
    assert(rbdeath > 0, 'This will affect B-Y retirement value if not true');

omega = 0.017;              %mortgage borrowing wedge
rmortgageAll = rbAll + omega;

cc_wedge = 0.103;
    rb_avg_wedge = zeros(Nb,1);
    rb_avg_wedge(1:bzeroInd-1) = cc_wedge;
    rb_avg_wedge_all = repmat(rb_avg_wedge,1,Na,Nz);

% 10Y TIPS Calibration
rate_r_OU = 0.2861;
sig_r_OU = 0.6280;

%Construct transition matrix (imposing reflecting barrier)
dr = 100*(rbAll(2)-rbAll(1));
la_mat_r = [100*(rbAll(1)-rmean)*rate_r_OU/dr - (sig_r_OU^2)/(dr^2) + sig_r_OU^2/(2*dr^2), -100*(rbAll(1)-rmean)*rate_r_OU/dr + sig_r_OU^2/(2*dr^2), 0, 0; ...
            sig_r_OU^2/(2*dr^2), 100*(rbAll(2)-rmean)*rate_r_OU/dr - sig_r_OU^2/(dr^2), -100*(rbAll(2)-rmean)*rate_r_OU/dr + sig_r_OU^2/(2*dr^2), 0;...
            0, 100*(rbAll(3)-rmean)*rate_r_OU/dr + sig_r_OU^2/(2*dr^2), -100*(rbAll(3)-rmean)*rate_r_OU/dr - sig_r_OU^2/(dr^2), sig_r_OU^2/(2*dr^2);...
            0, 0, 100*(rbAll(4)-rmean)*rate_r_OU/dr + sig_r_OU^2/(2*dr^2), -100*(rbAll(4)-rmean)*rate_r_OU/dr - sig_r_OU^2/(dr^2) + sig_r_OU^2/(2*dr^2)];
        
assert(sum(sum(la_mat_r,2)) < 1e-14);
assert(length(la_mat_r) == Nr);
assert(rbAll(2) - rbAll(1) == 0.01, 'Need to edit step-size in transition matrix');


%--------------------------------------------------------------------------
%% Blanchard-Yaari Retirement and Forced Adjustment
la_BY = 1/30;
la_forced = 1/15;


%--------------------------------------------------------------------------
%% HJB Equation Iteration Set-Up:
Delta     = 25; % step size in the implicit method
maxit     = 300;
tol       = 5e-5;
rng(1);
UseNoAdjustmentAsGuess = 1;
BilinW = 1; % when all neighboring points of (b',a') are non-adjustment points, 
            % use bilinear weights in dividing mass to be rebalanced across them in KF algorithm


%--------------------------------------------------------------------------
%% Create mesh: b in 1st dimension, a in 2nd, income in 3rd
[bbb, aaa, zzz] = ndgrid(b, a, zstates);
    aa_nob = squeeze(aaa(1, :, :));
    zz_nob = squeeze(zzz(1, :, :));
    bbb_negative = bbb.*(bbb<0);

M = Nb*Na*Nz;


%--------------------------------------------------------------------------
%% Build Exogenous Transition Matrices

% Matrix ZZ summarizing evolution of income
ZZ = [speye(Nb*Na)*la_mat_z(1,1), speye(Nb*Na)*la_mat_z(1,2), speye(Nb*Na)*la_mat_z(1,3);
    speye(Nb*Na)*la_mat_z(2,1), speye(Nb*Na)*la_mat_z(2,2), speye(Nb*Na)*la_mat_z(2,3);
    speye(Nb*Na)*la_mat_z(3,1), speye(Nb*Na)*la_mat_z(3,2), speye(Nb*Na)*la_mat_z(3,3)];


% Matrix AA summarizing evolution of a (no adjustment)
chi = -min(-xi*aaa,0)/da;
yy  = -max(-xi*aaa,0)/da + min(-xi*aaa,0)/da;
zeta= max(-xi*aaa,0)/da;

AA = sparse([]);
for nz = 1:Nz
    AAupdiag=zeros(Nb,1);
    for j=1:Na
        AAupdiag=[AAupdiag;zeta(:,j,nz)];
    end
    AAcentdiag= yy(:,1,nz);
    for j=2:Na-1
        AAcentdiag=[AAcentdiag;yy(:,j,nz)];
    end
    AAcentdiag=[AAcentdiag;yy(:,Na,nz)];
    AAlowdiag=chi(:,2,nz);
    for j=3:Na
        AAlowdiag=[AAlowdiag;chi(:,j,nz)];
    end
    AA = [AA; sparse(Nb*Na, Nb*Na*(nz-1)), ...
        spdiags(AAcentdiag,0,Nb*Na,Nb*Na)+spdiags(AAlowdiag,-Nb,Nb*Na,Nb*Na)+spdiags(AAupdiag,Nb,Nb*Na,Nb*Na), ...
        sparse(Nb*Na, Nb*Na*(Nz-nz))];
end


